Path integral virial estimator based on the scaling of fluctuation coordinates: 
Application to quantum clusters with fourth-order propagators 
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We first show that a simple scaling of fluctuation coordinates defined in terms of a given reference 
point gives the conventional virial estimator in discretized path integral, where different choices 
of the reference point lead to different forms of the estimator (e.g., centroid virial). The merit of 
this procedure is that it allows a finite difference evaluation of the virial estimator with respect 
to temperature, which totally avoids the need of higher-order potential derivatives. We apply this 
procedure to energy and heat capacity calculation of the (112)22 and Nei3 clusters at low temperature 
using the fourth-order Takahashi-Imada and Suzuki propagators. This type of calculation requires 
up to third-order potential derivatives if analytical virial estimators are used, but in practice only 
first-order derivatives suffice by virtue of the finite difference scheme above. From the application 
to quantum clusters, we find that the fourth-order propagators do improve upon the primitive 
approximation, and that the choice of the reference point plays a vital role in reducing the variance 
of the virial estimator. 



I. INTRODUCTION 



Imaginary time path integral provides a robust way for studying quantum statistical mechanics of many-particle 
systemsi^ In the framework of discretized path integral, this method maps a quantum system into multiple copies of 
virtual classical systems (called "beads") connected via harmonic springs. This isomorphism allows one to calculate 
structural and thermodynamic properties using conventional Monte Carlo or molecular dynamics methods. In practice, 
however, such calculation often becomes much more demanding than the classical counterpart, and thus a number 
of efficient techniques have been developed, e.g., collective sampling of multiple heads^*^ statistical estimators 
with low va ri anC eriL2i2ii2iiiii2ii2ii2ii&ii& and accurate approximations to the exact short-time propagator (or high- 
Oh; temperature density matrix), 3 ! 1 !! 18 , 13 ! 2 ", 21 ! 22 , 2 ^ 24 ! 25 ! 26 , 2 ?!^, 23 , 3 "!^, 32 !^!^, 3 ^ 

Our interest in this paper is in the latter two issues, namely the use of better statistical estimators and approxi- 
' mate propagators. Regarding the estimator, the most conventional path-integral estimators for internal energy are 
, thermodynamic^ and virial 7 " estimators. The former is obtained via direct temperature differentiation of the partition 
Q\ • function, while the latter is obtained by eliminating ill-behaved terms in the former through integration by parts. 
The virial estimator has an advantage that its variance is only weakly dependent on the number of beads, P, in 
contrast to the thermodynamic estimator whose variance grows linearly with P. We emphasize, however, that this 
reduction in the variance is achieved at the expense of using first-order potential derivatives that are absent in the 
. thermodynamic estimator. Although the first-order derivatives are often not a major computational problem, things 
become worse when one constructs a similar double virial estimator for heat capacity because it requires second-order 
potential derivatives. Despite this difficulty, the double virial estimator was used in a heat capacity calculation of 
water because other estimators exhibited too large statistical errors and could not be converged within simulation 
• \ time^£ To remedy this problem, Glaesemann and FriediS^ proposed a free-particle projection technique to reduce 
the variance of the thermodynamic estimator without using potential derivatives, and applied it to Arg clusters with 
considerable success at higher temperature. Predescu et alM> adopted a different strategy in their random series path 
O !■ integral (generalized form of the Fourier path integral), where they first scaled the amplitude of the Brownian bridge 
and then differentiated the scaled partition function via finite difference in order to obtain a viriallike estimator having 
no potential derivatives. With this method they calculated the quantum heat capacity of the Nei3 cluster at 4-14 K 
^ ■ with unprecedented accuracy^ 

Another issue that impacts the efficiency of path integral is the accuracy of approximate propagators. There exist 
a number of such approximations that aim at faster convergence to the P — > 00 limit than the standard primitive ap- 
proximation. In particular, the pair-product approximation 2,3 and the higher-order composite factorizationsiiiSiiSiS^ 
have proven to be successful in condensed-phase applications (see Ref. |27] for their useful comparisons) . The pair- 
product approximation replaces the exact high-temperature density matrix by the product of effective pairwise ones, 
and it has been shown to drastically reduce the number of beads for monoatomic fluidsi^* 3 ^ The fast convergence of 
this approximation was also exploited in semiclassical dynamical calculation of normal and superfluid helium. 37 While 
powerful for monoatomic fluids, the pair-product approximation becomes cumbersome when applied to molecular flu- 
ids because of the increased complexity of pair action. In this regard the higher-order propagators are appealing in 
that molecular fluids can be treated straightforwardly. In practice, however, the application of such propagators to 
molecular fluids is very scant compared to the primitive approximation. One reason may be that the higher-order 
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propagators involve the first-order potential derivatives, and the corresponding virial estimator for energy and heat 
capacity requires second- and third-order potential derivatives, respectively, resulting in a significant computational 
overhead. (Incidentally, Jang et al£& showed that for the Suzuki propagator the required order of potential derivatives 
can be reduced by using the virial theorem in operator form.) 

In this paper we present a method for evaluating the virial and double virial estimators in discretized path integral 
without using higher-order potential derivatives. This method is based on the coordinate scaling idea of Janke and 
SauesiH and the finite difference method of Predescu et al.i^ Specifically, we first show that a simple scaling of 
fluctuation coordinates defined in terms of a given reference point gives the conventional virial estimator, where 
different choices of the reference point lead to different forms of the estimator (e.g., centroid virial). This procedure 
reverts to the original coordinate scaling by Janke et al. when the reference point is set to the coordinate origin. We 
then take the temperature derivative of the scaled partition function by finite difference in order to avoid potential 
derivatives. We illustrate the above method by calculating energy and heat capacity of the (112)22 and Nei3 clusters at 
low temperature using the fourth-order composite propagators. This calculation requires up to third-order potential 
derivatives if analytical virial estimators are used, but in practice only up to first-order derivatives suffice by virtue of 
the finite difference scheme above. From the results of the application, we find that the fourth-order propagators do 
improve upon the primitive approximation, and that the choice of the reference point has a crucial role in reducing 
the variance of the virial estimator. 

The remainder of this paper is as follows: In Sec.[n]we describe the coordinate scaling and finite difference procedures 
mentioned above. In Sec. IIIII we apply the present method to the (H 2 )22 cluster at 6 K and Ne 13 cluster at 4-14 K 
and calculate their total energy, heat capacity, and distance distribution functions. Systematic comparisons are made 
among different types of propagators and estimators. In Sec. UVI we conclude. 



II. PATH INTEGRAL ESTIMATORS FOR ENERGY AND HEAT CAPACITY 



Conventional estimators 



We first summarize the conventional thermodynamic^ and virial^iS estimators for subsequent discussion. We 
suppose an /-dimensional system having the Hamiltonian H = T + V = 5Zf=i Pi I "^ m + ^( x ) with x = (xi, 
Using the primitive approximation to the canonical density operator, 



eV V 2 e - eT e- ey / 2 +0( £ 3 ), 
the partition function at inverse temperature (3 — 1/fcsT can be written as 



x f)- 
(1) 



Z((3) = tT{e- 0H ) = / dxi ■ • • / dx Pj0 (xi, . . . , x P ; 0) + 0{l/P z 



with 



p(xi,...,x P ;/3) 



mP 
2nh 2 (3 



Pf/2 



exp 



(2) 



(3) 



where x s is the system coordinate in the s-th time slice (or "beads") with xo = xp. The thermodynamic estimator 
is obtained by direct temperature differentiation of Eq. (J2J : 



E(f3) 



1 dz(p) 

Z(f3) 8f3 
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where (• ■ •) denotes an ensemble average over the sampling function p(x.\, . . . , xp; (3). The drawback of this estimator 
is that its variance grows with P due to cancellation of the first two terms in the right-hand side of Eq. ©. This 
difficulty can be avoided by using the relation, 
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,xp;/3) = ~(P-g)f J dxi ■ ■• J dx P p(x 1 , . . . ,x P ;/3), (6) 
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which arises from integration by parts. In Eq. JJJ, x* is a given "reference" point and g is a constant that depends 
on the definition of x*. In this paper we consider three choices of x*, namely x* = 0, xp, and x c , where x c is the 
centroid of the imaginary-time path given by 

I p 

s=l 

With these choices the value of g becomes^ 

?' X I = °' j (8) 
1, x = Xp and x c . v ' 

Because the kinetic action in Eq. J3J is doubled by the "virial operator" in the square bracket in Eq. the 
following path integral virial theorem holds: 

p p 

mP 



\ " S — 1 S— 1 / 

Eliminating the first two terms in Eq. JSJ through the above relation, we have the following virial estimator for energy: 



1 \- 



fg 

2/3 ' P 



s=l 



(10) 



For convenience we will refer to the above estimator with x* = 0, xp, and x c as the origin-, bead-, and centroid- 
reference virial estimators, respectively. We note that the origin-reference virial estimator gives an incorrect result for 
unbounded systems^ (e.g., ey vanishes for a free particle) although the bead- and centroid-reference virial estimators 
remain valid. The reason is that in the former the integral of p(xi, . . . , xp; (3) over the whole coordinate space 
is divergent, which invalidates Eq. JfjJ, while in the latter the integration by parts can be performed in one less 
dimensions with some coordinate fixed (e.g., xp in the bead-reference virial). Despite this deficiency, the origin- 
reference virial estimator can be applied to quantum clusters if the contribution of the center of mass is properly 
taken into account piLiii 

Heat capacity estimators can be obtained in a similar manner and are fully described in Ref. I13L The resulting 
double thermodynamic estimator contains no potential derivatives but its variance grows rapidly as P 2 . The double 
virial estimator has a favorable variance weakly dependent on P but it requires second-order potential derivatives, 
resulting in an increased computational effort^ 



B. Virial estimator via coordinate scaling 

The virial estimator in Eq. (|10(l can also be obtained by a scaling of fluctuation coordinates as mentioned in the 
Introduction. This is achieved by first considering the partition function at a different temperature /3': 

Z((]') = fdx; ••• J dx P p(xi, . . . ,x P ;/3'), (11) 

where p is the density function in Eq. @- To eliminate ill-behaved terms in the thermodynamic estimator, we 
introduce a new set of variables (xi, . . . , xp) as 

'|(x, -x). (12) 

for s = 1, ... , P. x* in Eq. (|12|) is a reference point that has the same meaning as in the preceding section, i.e., x* = 0, 
xp, or x c . The Jacobian of this transformation is 

/0t\(P-9)f/2 

dx.^ ■ ■ ■ dx.' P = ( — I g?xi ■ • • dxp, (13) 
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where g is given by Eq. JSJl. Since the transformation in Eq. I|12[) suggests 

P s=l ^ s=l 

the partition function in Eq. may be written as 

Z{0) = J dxx--- J dxpp( Xl ,...,xp;^W) 



with 



r(p') = (|) /9/2 ex p |-p X> v ( x *) - /^)]} 



Using the above equation the internal energy is obtained as follows, 



E(I3) = - 
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(18) 



0'=0 



where we have explicitly denoted the /3'-dependence of x' s . ev in Eq. (|18(1 can be shown identical to that in Eq. I jlOjl 
by taking the /^'-derivative analytically. Instead, we may take the /3'-derivative via finite difference in order to avoid 
potential derivatives*^ 



V 2ft 2P5ft 



i p 
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Similarly, the constant volume heat capacity, 

dE{T) 2 
Cv(ft) = dT = k B ft 

can be obtained using the following expression, 
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(19) 

(20) 
(21) 
(22) 



E(ft) 



(23) 



where we may use finite difference to evaluate the second derivative with respect to ft' . Although there are other 
schemes for performing finite difference, e.g., 

' R(ft + 8ft) - R(ft - 8ft) \ 
28ft /' 

it has a narrower range of acceptable values of 8 ft than Eq. (|19fl due to the exponential behavior of R(ft'). Therefore 
we will use only Eqs. (|19[) and i122[) in the following sections. 

We emphasize that the finite difference scheme above is qualitatively different, e.g., from that performed for internal 
energy with respect to temperature, 



Cv(ft) 



E(T + ST) - E(T - 8T) 
28T ' 



(24) 



where E(T) contains statistical error and ST must be taken sufficiently large so that \E(T + ST) — E(T — ST) \ is much 
larger than the statistical error in E(T ± ST). On the other hand, the finite difference in Eq. I|19fl is performed for a 
statistical error-free quantity, ft'V(x' s (ft')), and thus 8 ft can be taken as small as machine precision allows. In practice, 
however, acceptable values of 5ft may depend on the stiffness of the potential as well as thermodynamic conditions 
under study (e.g., particle density), so one needs to check the convergence by repeating a very short simulation with 
different values of 8ft. Our typical choice of 8ft is 10 _4 /3 (see Sec. Illlll . 
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C. Using fourth-order composite propagators 



An appealing feature of the finite-difference scheme in Sec. Ill Bl is that it does not require potential derivatives 
higher than those existing in the discretized action. This means that energy and heat capacity can be calculated 
with no potential derivatives when the primitive approximation is used, and only up to first-order derivatives are 
needed when the fourth-order composite propagators are used. The generalized SuzukiiLiS^ and Takahashi-Imadai^ 
approximations fall into the latter category. The Suzuki approximation factorizes the exact short-time propagator as 

e . 2eH = e _ £ y e/ 3 e - £ T e -4 £ y m /3 e - e T e - e y c /3 + (25) 

where V m and V e are effective potentials that involve first-order potential derivatives (see Refs. and|2j| for details). 
With this factorization the approximate partition function becomes 

Z((3) = y"dx 1 .-.y"^P j0 (4) (x 1 ,...,xp;/3)+O(l/P 4 ) (26) 

with 

p W { ^...^P)=(j^ Pi '\^ (27) 

where V a is a time-slice dependent effective potential defined by 

fc(x; (3) = V(x) + 4(/3/P) 2 C(x) (28) 

with 

dV(x) 



C(x) = [V,[T,V\] = - 
m 



(29) 



while w s and d s are a set of coefficients given by 

2/3, s = even 



4/3, s = odd, (30) 



and 



J a/6, s = even, , , 

fls - \ (l-a)/12, s = odd, 

where a is an arbitrary parameter within [0,1]. The partition function for the Takahashi-Imada approximation^ can 
also be expressed in the form (|27|) with w s = 1 and d s = 1/24, although this approximation is not based on a genuine 
factorization such as Eq. The fourth-order approximation in Eq. Il26[l differs from the primitive, second-order 

one only in that the bare potential is replaced by the slice-dependent effective potential in Eq. I|28|l. and that the 
weight factors {w s } are introduced in a way similar to Simpson's quadrature rule, ft is thus straightforward to apply 
the procedure in Sec. Ill Bl to obtain finite-difference virial estimators having no higher-order potential derivatives. 
Specifically, the statistical average is now taken over in Eq. J57J), and we modify R{(3') in Eq. I|16|) as 

R W)=(j) fa ' e X p|-i^^[/3'y s (x^;/3')-/3K(x s ^)]|, (32) 
and the virial estimators in Eqs. (|18(l and Ill'2|) as follows: 
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TABLE I: Internal energy (K/molecule) calculated using the centroid virial estimator based on the finite difference scheme in 
Sec. Ill Bl PA, TIA, and SA denote the primitive, Takahashi-Imada, and Suzuki approximations, respectively, a is an arbitrary 
parameter within [0,1] involved in the Suzuki approximation. The figure in parentheses is one standard deviation on the last 



digit. 


P 


PA 


TIA 


SA(q = 0) 


SA (q=1/2) 


SA(a=l) 


20 


-27.52(1) 


-21.80(1) 


-23.24(1) 


-20.94(1) 


-21.66(1) 


40 


-21.54(1) 


-18.80(1) 


-19.41(1) 


-18.34(1) 


-18.38(1) 


60 


-19.78(1) 


-18.13(1) 


-18.41(1) 


-17.86(1) 


-17.78(1) 


80 


-18.98(1) 


-17.90(1) 


-18.08(1) 


-17.72(1) 


-17.61(1) 


100 


-18.56(1) 


-17.79(1) 


-17.90(1) 


-17.66(1) 


-17.60(1) 


120 


-18.34(1) 


-17.76(1) 


-17.84(1) 


-17.68(1) 


-17.63(1) 


160 


-18.09(1) 


-17.73(1) 


-17.76(1) 


-17.68(1) 


-17.65(1) 



TABLE II: Heat capacity (in unit of fcs) calculated using the double centroid virial estimator based on the finite difference 
scheme in Sec. Ill Bl Other details are the same as in Table Q 



p 


PA 


TIA 


SA(a = 0) 


SA (a=l/2) 


SA(a=l) 


20 


80.6(4) 


59.3(5) 


65.5(5) 


58.8(5) 


64.2(5) 


40 


55.5(4) 


44.5(5) 


47.4(4) 


42.5(5) 


44.4(5) 


60 


47.7(4) 


38.9(4) 


40.9(4) 


36.9(5) 


36.4(5) 


80 


42.5(4) 


37.7(4) 


38.5(4) 


35.8(4) 


35.1(5) 


100 


40.7(4) 


35.8(4) 


37.3(4) 


35.2(4) 


34.5(4) 


120 


39.3(4) 


35.7(4) 


35.6(4) 


34.9(4) 


34.8(4) 


160 


37.6(4) 


34.6(4) 


35.5(4) 


34.3(4) 


33.6(4) 



III. APPLICATION TO QUANTUM CLUSTERS 
A. (H 2 )22 cluster at 6 K 

We illustrate the above procedure by first calculating the energy and heat capacity of the (1-12)22 cluster at 6 K. 
The physical model is identical to that used in the previous studies ^^iJ^ Briefly, the system potential consists 
of Lennard- Jones (LJ) pair interactions with £lj = 34.2 K and ctlj = 2.96 A, where the hydrogen molecules are 
treated as distinguishable spherical particles with their mass being 2 amu. Since a cluster in vacuum at any positive 
temperature is metastable with respect to evaporation, a confining potential of the form, 

h_JM , ,35, 

i— 1 x ' 

is added to the sum of LJ potentials to prevent any molecules from permanently leaving the cluster. In Eq. (|35|l . is 
the position of particle i, R is the center of mass of the cluster, TV is the number of particles, and R c is the confining 
radius chosen as 4ctlj- 

Statistical sampling of imaginary-time paths was performed with Monte Carlo (MC) methods, where one cycle is 
defined such that each particle is moved once on average by the staging algorithm. The implementation is the 
same as described in Ref. 0- The staging length j (the number of beads that are collectively moved) was determined 
by adjusting the acceptance ratio to 50 %, which resulted in j w P/A regardless of the value of P. In addition to the 
staging move, we applied the whole-chain move^ every two cycles to accelerate the statistical convergence. A single 
run consisted of 5 x 10 5 cycles for equilibration followed by 4 x 10 6 cycles for data accumulation, which took several 
days for P = 160 using a Pentium4 3.8 GHz PC. 

Tables [Q and CU list the energy and heat capacity obtained using the primitive, Takahashi-Imada, and Suzuki 
approximations. Figure ^ illustrates the systematic convergence of those values to the P — > 00 limit. The centroid- 
reference virial estimator was used throughout based on the finite difference scheme presented in Sec. Ill Bl The 
relative statistical error was estimated to be on the order of 0.1 and 1 % for energy and heat capacity, respectively, by 
using a blocking procedure with 2000 blocks each of 2000 cycles. The present calculation needed only up to first-order 
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FIG. 1: Systematic convergence of (a) energy (K/molecule) and (b) heat capacity (in unit of fcs) of the (112)22 cluster at 6 K 
as a function of the Trotter number P. 




FIG. 2: Statistical error in (a) energy (K/molecule) and (b) heat capacity (in unit of ks) of the (112)22 cluster at 6 K as a 
function of the Trotter number P. Five different estimators are compared: the thermodynamic estimator, the virial estimator 
with different choices of the reference point, and a modified virial estimator in Eq. Ij36ft . Errors in the heat capacity were 
estimated using the prescription given in Ref. HH 



potential derivatives as mentioned in the Introduction. Acceptable values of the stepsize 8(3 ranged broadly from 
10 _3 /3 to 10 _6 /3, where the smallest value was determined by round-off errors in the heat capacity. In this paper we 
set 8/3 to 10~ 4 /9, which practically gave the same result as when the analytical virial estimator was used. 

Figure^shows that the fourth-order approximations improve remarkably upon the primitive approximation for both 
the energy and heat capacity. For example, to achieve a systematic error in the energy less than 0.25 K/molecule j^^i^ 
the primitive approximation requires P > 200 while the Suzuki approximation having a > 0.5 attains the same 
accuracy with P = 60, thus reducing the necessary value of P by a factor of ~3. This acceleration of systematic 
convergence is similar to that observed by Brualla et id^L in the study of liquid 4 He at 5.1 K using the Takahashi- 
Imada approximation. Regarding the converged values of internal energy, the present result (E = —17.68 ± 0.01 
K/molecule) obtained using the Suzuki approximation with P = 160 and a = 0.5 is in excellent agreement with 
the most accurate estimate (E = —17.69 ± 0.01 K/molecule) obtained by Predescu et alSk using the Wiener-Fourier 
reweighted path integral method. Comparing different fourth-order propagators, we see that the Suzuki propagator 
with a = 0.5 and 1.0 converges somewhat faster than that with a — 0.0 or the Takahashi-Imada approximation. 

Figure [2] compares statistical errors in the energy and heat capacity obtained with different estimators. The 
discretization was performed using the Suzuki approximation with a — 0.5. This figure shows that the thermodynamic 
estimator has a growing variance with P while the virial estimators have a nearly constant variance. We should note, 
however, that the variance of the virial estimator strongly depends on the choice of the reference point x* in Eq. i|12[l . 
That is, the origin-reference virial estimator exhibits a significantly larger variance than the bead- or centroid-reference 
estimators, indicating that it is more advantageous to choose x* in Eq. I|12l) as xp or x c than the coordinate origin. 
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FIG. 3: (a) Pair radial distribution function and (b) distance distribution function from the center of mass of the (¥[2)22 cluster 
at 6 K. SA and PA denote the quantum results obtained with the Suzuki and primitive approximations, respectively. The 
classical result is plotted in short dashed line (with its height scaled by a factor of 1/2 to fit in the panel). 



Figure also plots the statistical error of the modified centroid virial estimator that may be used in path integral 
molecular dynamics In the latter method the fourth-order composite propagators become expensive if in 
Eq. I|27l) is used directly as a sampling function, because the "forces" exerted on the beads require second-order 
potential derivatives. This problem can be avoided, for example, by excluding the force square terms in Eq. (|29[) from 
the sampling function^ The resulting modified expression for the energy is 

with 

(37) 



Ap = cxp I - w * d s (P/P) S C(x s ) I 



where various symbols are the same as in Sec. Ill (Jl (• • -) mo d in Eq. (|36|) denotes an ensemble average over the following 
sampling function: 

Pmod(xi,...,x P ) = expi -^r^J2(*s -Xs-i) 2 - ^J2w s V(x s ) 1 . (38) 

L ^ s—l s—1 ) 

This method gives the true expected value of energy as P is increased, but the statistical error becomes larger than 
the original scheme in Sec.|n] 23 Figure [21 shows that the variance obtained with this method is quite large for small 
values of P but is reduced to a manageable size if P is increased to >40. Thus, excluding force square terms from the 
sampling function seems a viable option if molecular dynamics methods are used as a statistical sampler. 

Figure |3] illustrates the classical and quantum results of the pair radial distribution function and the distance 
distribution function from the cluster center of mass defined by 

PpairMcx/f^^r-rgM, (39) 

\ i<j s I 

and 

Pcm (r)oc(^^«5(r-ArW)y (40) 



respectively, where = | r i — r j I and Ar^ = |r s ; — R| . The sum over time slices is performed for all (only even) values of 
s when the primitive (Suzuki) approximation is used^ Here we do not consider the Takahashi-Imada approximation 
because it requires a nontrivial modification to the estimator^ We see from Fig. [3] that the classical cluster has 
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a rigid, solidlike structure at this temperature^ while the quantum cluster has a liquidlike structure due to large 
zero-point energies and tunneling effects. This figure also shows that both the primitive and Suzuki approximations 
with P = 20 already give a good approximation to the practically exact result obtained with P = 160, indicating that 
structural properties converge much faster than the energy and heat capacity as a function of P. 



Nei3 is one of the smallest clusters that exhibit solid-liquid-like (or melting) transition, and it has been studied 
extensively using a variety of theoretical methods 4^i^ 4 ^^S^£' 4 ^ The classical melting point is located at around 
10 K and it is lowered by about 10 % due to prominent quantum effects. The heat capacity is a useful quantity for 
characterizing such a cluster phase transition. Neirotti et alw^ calculated the heat capacity of Nei3 using a double 
virial estimator (in analytical form) designed for the Fourier path integral, but unfortunately their results exhibited 
a large statistical error (about 10 fee) in the low temperature region. Predescu et alw^ calculated the same quantity 
using their viriallike estimator (in finite-difference form) in the framework of random series path integral, and as 
mentioned in the Introduction they obtained highly converged results with statistical errors less than 1 ks- Because 
the two calculations used the same number of Monte Carlo samples, this reduction in statistical error corresponds 
roughly to 100 times acceleration in convergence rate. Then a natural question that arises is what is the dominant 
factor that reduced the statistical error. We find, however, that this question is rather difficult to answer because 
there are quite a few technical differences in their calculations. 

As such, to get some insights into the above question, we have re-calculated the heat capacity of Nei3 using the 
discretized path integral. The computational details are basically the same as in the preceding section, and the relevant 
parameters were set as closely as possible to those in Refs. Ill andlTEl Specifically, the Lennard- Jones parameters were 
set to £lj = 35.6 K and c7lj = 2.749 A, and the mass of Ne was 20.0 amu. The confining radius R c in Eq. I|35fl was 
chosen as 2olj- The number of Monte Carlo cycles was 4 x 10 6 . We also performed the replica-exchange (or parallel 
tempering) Monte Carloi2iii2i£i to avoid nonergodicity problem at low temperature. The number of replicas was set 
to 21, and the replica temperatures were distributed over the interval [4,14] K with even spacing. The exchange move 
was attempted every 10 Monte Carlo cycles. This setting ensured the acceptance ratio of exchange moves to be >10 



Figure^lplots the heat capacity thus obtained as a function of temperature. Four combinations of the approximate 
propagator and the heat capacity estimator are examined, namely: 

(a) primitive approximation + double thermodynamic estimator; 

(b) primitive approximation + origin-reference double virial estimator; 

(c) primitive approximation + bead-reference double virial estimator; 

(d) Suzuki approximation (a = 0.5) + centroid-reference double virial estimator. 

In cases (b), (c), and (d) the virial estimator was evaluated using the finite difference scheme in Sec. Ill Bl while in 
case (a) the heat capacity was calculated using the double thermodynamic estimator: 



where Et is given by Eq. (JSJ. In all cases the number of time slices was set to P = 8-24. Also shown in Fig. 01 is 
the highly accurate results obtained by Predescu et alr& (circles) and the classical heat capacity (dotted line). This 
figure reveals that the origin-reference virial estimator has much larger statistical errors than the other three cases. 
Also interesting is the fact that the variance of the double thermodynamic estimator is rather small and close to the 
bead- and centroid-reference double virial estimators. This tendency is qualitatively similar to that observed for the 
heat capacity of the hydrogen cluster in Fig. [21 where the origin-reference virial estimator has the largest statistical 
error in the small P region. Regarding the systematic convergence to the P — > oo limit, Fig. 0] (d) shows that the 
Suzuki propagator again provides a noticeable improvement over the primitive approximation, and that P = 20 is 
sufficient to reach systematic convergence within 1 ks- 

What is more important about the question discussed above is that cases (b) and (c) correspond qualitatively to 
the calculation by Neirotti et alm^ and Predescu et al.^ respectively. More precisely, the double virial estimator of 
Neirotti et al. may be regarded as origin-reference because their estimator vanishes when the interaction potential is 
set to (and thus the contribution of the center of mass was treated separately), while the finite-difference estimator 
of Predescu et al. may be viewed as bead-reference because the relevant Brownian bridge is defined in terms of a set 



B. Neis cluster at 4-14 K 
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FIG. 4: Heat capacity of the Nei3 cluster as a function of temperature. SA and PA denote the quantum results obtained with 
the Suzuki and primitive approximations, respectively. Different estimators are used in each panel: (a) double thermodynamic 
estimator in Eq. 14H : (b) origin-reference double virial estimator; (c) bead-reference double virial estimator; (d) centroid- 
reference double virial estimator. The statistical errors in panels (a), (c), and (d) are comparable to the width of the line while 
that in panel (b) is about 5 ks in the low temperature region. The highly accurate results obtained by Predescu et al. fRef. ll5l) 
are plotted by circles. The classical result is plotted in dotted line. 



of "physical coordinates" (equivalent to a single bead). Thus, we think that the dominant factor that made a large 
difference in their calculations is the choice of the reference point in the virial estimator, rather than whether the 
estimator was evaluated analytically^ or numerically via finite difference^ if we consider the fact that cases (b) , (c) , 
and (d) above were treated using the finite difference method in Sec. Ill Bl 

Finally, Fig. El illustrates the pair distribution functions at 4 and 10 K, which are very similar to those presented 
in Refs.|4a andl47l The cluster takes a solidlike structure at 4 K while it starts to form a liquidlike structure at 10 
K (slightly above the melting temperature). Comparing the classical and quantum results, we see that the positions 
of the classical peaks are shifted outward and their widths broadened when quantum effects are made operative. The 
degree of broadening is much smaller than that observed for the hydrogen cluster in Fig.[3]due to quasiclassical nature 
of Nei3. Figure[S]also shows that the primitive and Suzuki approximations with P — 24 give almost indistinguishable 
results, verifying the fast convergence of structural properties with respect to P. 



IV. CONCLUSIONS 



In this paper we have presented a coordinate scaling procedure for obtaining the conventional virial estimator and 
discussed its efficient evaluation using finite difference with respect to temperature. This procedure allowed us to 
apply the fourth-order propagators to quantum clusters using only the first-order potential derivatives. From the 
results of the application, we find that setting the reference point in the virial estimator to the coordinate origin (the 
path centroid) gives the largest (smallest) statistical errors. This result is also in qualitative agreement with previous 
studies on Ar clusters and liquid water >i2iiii& Thus, despite its extensive use in the literature, it is not recommended 
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FIG. 5: Pair radial distribution functions of the Nei3 cluster at (a) 4 K and (b) 10 K. SA and PA denote the quantum results 
obtained with the Suzuki and primitive approximations, respectively. The classical result is plotted in dashed line. 



to use the origin-reference virial estimator in quantum clusters and condensed phase systems because of the large 
variance as well as unnecessary complication due to unbounded degrees of freedom. 

We end this paper by mentioning some possible application of the present method. One example is a short-time 
approximation to the quantum correlation function Cab^), e.g., 

Cab® = C AB (0) + ^C AB (Q)t 2 + ■■■ ~ C AB (0)exp 

Taking the real-time derivatives of Cab (t) at t = along the imaginary-time axis results in a path integral calcula- 
tion similar to that of heat capacity,? 3 which implies a similar reduction in statistical errors via coordinate scaling. 
Such an idea is particularly relevant, e.g., to an approximate calculation of chemical reaction rates^ 3 ^*^*^ or vibra- 
tional relaxation rates^i Another possible application is the path integral ground state (or variational path integral) 
methodsjSiSiS where several different estimators arise in natural analogy to finite temperature path integral. 



1 C AB (0) 2 
2C AB (0) 



(42) 
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